source("code/invertRIconfInt.R")
source("code/seqblock1.R")
source("code/seqblock2k.R")
source("code/mahal.R")

load("data/SMRPreplication.RData")

smrp.no <- x[x$withdrawn == "no", ]
smrp.no$tr <- as.numeric(smrp.no$condition)-1

## Create inverted randomization inference confidence interval (optional -- loaded below):
set.seed(742)
out <- invertRIconfInt(dat = smrp.no, outcome.var = "pssi2_totFake", tr.var = "tr", tau.abs.min = -6, tau.abs.max = 6, tau.length = 100, n.sb.p = 1000)
load("data/outRIfakeFULLfine.RData")

dm.ci.95 <- t.test(pssi2_totFake ~ condition, data = smrp.no, conf.level = .95)$conf.int
dm.ci.90 <- t.test(pssi2_totFake ~ condition, data = smrp.no, conf.level = .90)$conf.int
dm.ci.80 <- t.test(pssi2_totFake ~ condition, data = smrp.no, conf.level = .80)$conf.int

## Figure 9 (right):
ysRI <- rep(3/4, 3)
ystt <- rep(1/4, 3)
plot(-100,-100, xlim=c(-8,8), ylim=c(0,1), xlab="(Adjusted) Outcome Estimate", ylab="", axes = FALSE)
axis(1, at = c(-8, -4, 0, 4, 8))
axis(2, at = c(3/4, 1/4), labels = c("Rand Inf 
                                     conf ints ", "t-test   
                                     conf ints"), lty=0, las = 1, cex.axis=.9)
segments(sort(out$ci95)[2], ysRI[3], max(out$ci95), ysRI[3], lwd = 5, col="lightgrey")
segments(min(out$ci90), ysRI[2], max(out$ci90), ysRI[2], lwd=10, col="grey")
segments(min(out$ci80), ysRI[1], max(out$ci80), ysRI[1], lwd=20)
segments(dm.ci.95[1], ystt[3], dm.ci.95[2], ystt[3], lwd = 5, col = "lightgrey")
segments(dm.ci.90[1], ystt[2], dm.ci.90[2], ystt[2], lwd = 10, col = "grey")
segments(dm.ci.80[1], ystt[1], dm.ci.80[2], ystt[1], lwd = 20)

## Values:
summary(out$ci95)
summary(out$ci90)
summary(out$ci80)

## How much narrower?
## 9-15\% narrower:
((dm.ci.95[2]-dm.ci.95[1])- (max(out$ci95)-min(out$ci95)))/(dm.ci.95[2]-dm.ci.95[1])
((dm.ci.90[2]-dm.ci.90[1])- (max(out$ci90)-min(out$ci90)))/(dm.ci.90[2]-dm.ci.90[1])
((dm.ci.80[2]-dm.ci.80[1])- (max(out$ci80)-min(out$ci80)))/(dm.ci.80[2]-dm.ci.80[1])
